#ifdef _OPENMP 
#include <omp.h> 
#endif
#include <iostream>
#include<stdlib.h>
#include<stdio.h>
#include<time.h>
#include <arm_neon.h>
using namespace std;
const int N = 1000;
const int NUM_THREADS=7;
float m[N][N];
int i, j, k;
float tmp;
void m_reset() {
	for (int i = 0; i < N; i++)
	{
		for (int j = 0; j < i; j++)
			m[i][j] = 0;
		m[i][i] = 1.0;
		for (int j = i + 1; j < N; j++)
			m[i][j] = rand();
	}
	for (int k = 0; k < N; k++)
		for (int i = k + 1; i < N; i++)
			for (int j = 0; j < N; j++)
				m[i][j] += m[k][j];
}
void serial() {
	for (int k = 0; k < N; k++)
	{
		for (int j = k + 1; j < N; j++)
			m[k][j] = m[k][j] / m[k][k];
		m[k][k] = 1.0;
		for (int i = k + 1; i < N; i++)
		{
			for (int j = k + 1; j < N; j++)
				m[i][j] = m[i][j] - m[i][k] * m[k][j];
			m[i][k] = 0;
		}
	}

}
void omp_LU() {
#pragma omp parallel num_threads(NUM_THREADS), private(i, j, k, tmp) 
	for (k = 0; k < N; ++k) {
		// 串行部分，也可以尝试并行化
	#pragma omp single
		{
			tmp = m[k][k];
			for (j = k + 1; j < N; ++j) {
				m[k][j] = m[k][j] / tmp;
			}
			m[k][k] = 1.0;
		}
		// 并行部分，使用行划分
	#pragma omp for
		for (i = k + 1; i < N; ++i) {
			tmp = m[i][k];
			for (j = k + 1; j < N; ++j) {
				m[i][j] = m[i][j] - tmp * m[k][j];
			}
			m[i][k] = 0.0;
		}
	}
}
void omp_SIMD() {
#pragma omp parallel num_threads(NUM_THREADS), private(i, j, k, tmp) 
	for (k = 0; k < N; ++k) {
		
	#pragma omp single
		{
			tmp = m[k][k];
			float32x4_t vt=vld1q_f32(&tmp);
        for(int j=k+1;j<N-3;j+=4){
            float32x4_t va=vld1q_f32(&m[k][j]);
            va=vdivq_f32(va, vt);
            vst1q_f32(&m[k][j],va);
        }
        for(;j<N;j++)
            m[k][j]=m[k][j]/tmp;
		m[k][k] = 1.0;
		}
		// 并行部分，使用行划分
	#pragma omp for
    for(int i=k+1;i<N;i++){
		int a=0;
        float32x4_t vaik;
        vaik=vld1q_dup_f32(m[i]+k);
        tmp = m[i][k];
        for (int j = k + 1; j +4< N; j+=4,a=j)
        {
          
            float32x4_t vakj=vld1q_f32(m[k]+j);
            float32x4_t vaij=vld1q_f32(m[i]+j);
            float32x4_t vx=vmulq_f32(vakj,vaik);
            vaij= vsubq_f32(vaij,vx);
            vst1q_f32(m[i]+j,vaij);
        }
			 for(int j=a+1;j<N;j++)
				m[i][j] = m[i][j] - tmp * m[k][j];
			
			m[i][k] = 0.0;
		}
    }
}
    

int main() {
	
	m_reset();
    struct timespec sts,ets;
    timespec_get(&sts,TIME_UTC);
    serial();
    timespec_get(&ets,TIME_UTC);
    time_t dsec=ets.tv_sec-sts.tv_sec;
    long dnsec=ets.tv_nsec-sts.tv_nsec;
    if (dnsec<0){
             dsec--;
    dnsec+=1000000000ll;
    }
    cout<<"LU分解单线程:";
    printf ("%lld.%09llds\n",dsec,dnsec);
	// 在外循环之外创建线程，避免线程反复创建销毁，注意共享变量和私有变量的设置

	m_reset();
	timespec_get(&sts,TIME_UTC);
	omp_LU();
	timespec_get(&ets,TIME_UTC);
    dsec=ets.tv_sec-sts.tv_sec;
    dnsec=ets.tv_nsec-sts.tv_nsec;
    if (dnsec<0){
             dsec--;
    dnsec+=1000000000ll;
    }
	cout << "LU多线程:";
    printf ("%lld.%09llds\n",dsec,dnsec);

	m_reset();
	timespec_get(&sts,TIME_UTC);
	omp_SIMD();
	timespec_get(&ets,TIME_UTC);
    dsec=ets.tv_sec-sts.tv_sec;
    dnsec=ets.tv_nsec-sts.tv_nsec;
    if (dnsec<0){
             dsec--;
    dnsec+=1000000000ll;
    }
	cout << "SIMD多线程:";
    printf ("%lld.%09llds\n",dsec,dnsec);
}